% Simulating Crestor market share under different scenarios
% code skeleton based on smle_pairwise_generic_correlated
clear
clc


%% Take parameter estimates from smle_pairwise_generic_correlated
delta = zeros(4,1); % baseline is 0.25 split between 3 brands and outside option
delta(4) = 0.5;
gamma = 5 * eye(4,4);
gamma(4,1) = 5;


%% Play out market shares over some interval using 
% delta = -inf if the drug is off the market
shares_init = [0.4;0.4;0;0];
g = 0.02;

% simulate from 2002 to 2011
% y=1 is 2002Q1
Y = 40;
delta_reference = repmat(delta,1,Y);




%% Scenario 1: baseline
delta_matrix = delta_reference;
delta_matrix(3,1:6) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_basic = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_basic(:,y) = shares_curr;
end



%% Scenario 2: start of 2003 year earlier
delta_matrix = delta_reference;
delta_matrix(3,1:4) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_early = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_early(:,y) = shares_curr;
end


%% Scenario 3: enter 1 year later
delta_matrix = delta_reference;
delta_matrix(3,1:8) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_later = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_later(:,y) = shares_curr;
end


%% Scenario 4: enter 2 years later
delta_matrix = delta_reference;
delta_matrix(3,1:12) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_ll = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_ll(:,y) = shares_curr;
end

%% Scenario 5: enter at the same time as generic (2006)
delta_matrix = delta_reference;
delta_matrix(3,1:16) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_lll = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_lll(:,y) = shares_curr;
end


%% Scenario 5: enter after the generic (2007)
delta_matrix = delta_reference;
delta_matrix(3,1:20) = -Inf;
delta_matrix(4,1:18) = -Inf;

shares_llll = zeros(4, Y);
shares_curr = shares_init;

for y=1:Y
   curr_delta = delta_matrix(:,y);
   shares_curr = demand(curr_delta, gamma, shares_curr, g);
   shares_llll(:,y) = shares_curr;
end


h=figure;
hold on
t=2002:0.25:2011.9;

plot(t,shares_basic(1,:),'b',t,shares_basic(2,:),'b:',t,shares_basic(3,:),'b--', t,shares_basic(4,:),'b-.') % columns of Y vs. columns of X


set(gca,'FontSize',16)
set(gcf,'color','w');
xlabel('Year','FontSize',18)
ylabel('Market Share','FontSize',18)
legend({'Zocor','Lipitor','Crestor','Zocor (g)'},'FontSize',12)

saveas(h, 'sim_all_drugs.png');
savefig('sim_all_drugs.fig');
hold off




h=figure;
hold on
t=2002:0.25:2011.9;

plot(t,shares_basic(3,:),'b',t,shares_early(3,:),'b:',t,shares_later(3,:),'b--', t,shares_ll(3,:),'b--', t,shares_lll(3,:),'b--', t,shares_llll(3,:),'b--') % columns of Y vs. columns of X
%legend('Baseline','2002','2004','2005','2006','2007')

set(gca,'FontSize',16)
set(gcf,'color','w');
xlabel('Year','FontSize',18)
ylabel('Market Share','FontSize',18)

saveas(h, 'crestor_sim.png');
savefig('crestor.fig');
hold off






